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Abstract 

This paper is concerned with stochastic reaction-diffusion kinetics governed by the 
reaction-diffusion master equation. Specifically, the primary goal of this paper is to 
provide a mechanistic basis of Turing pattern formation that is induced by intrinsic 
noise. To this end, we first derive an approximate reaction-diffusion system by us- 
ing linear noise approximation. We show that the approximated system has a certain 
structure that is associated with a coupled dynamic multi-agent system. This obser- 
vation then helps us derive an efficient computation tool to examine the spatial power 
spectrum of the intrinsic noise. We numerically demonstrate that the result is quite 
effective to analyze noise-induced Turing pattern. Finally, we illustrate the theoretical 
mechanism behind the noise-induced pattern formation with a I-L2 norm interpretation 
of the multi-agent system. 

1 Introduction 

The auto-regulation mechanism of the biological pattern formation is a long standing ques- 
tion in biology. It was Turing ^ who first presented a mathematical model that possibly 
accounts for the pattern formation in embryonic development. Specifically, it was shown 
that a pair of reaction-diffusion equations can spontaneously exhibit spatial structure after 
a small perturbation to a spatially homogeneous equilibrium. Nowadays, such patterns are 
called Turing pattern, and analysis ranges to a wide variety of biological applications (see 
[2] for example). 

Although many existing studies rely on deterministic reaction-diffusion models, intrinsic 
noise is often a unignorable factor leading to a drastic change of the dynamics [3 . It is 
known that the dynamics of stochastic chemical reactions in a cell follows the chemical 
master equation (CME) [4 . A standing assumption of the CME is that reactants are well- 
mixed in a cell. However, recent studies revealed that molecules can localize inside a cell, 
and the localization plays an important role, for example, in cell division [5 . Hence, it is 
important to analyze the stochastic dynamics of biochemical reactions under the spatially 
inhomogeneous environment. 
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The reaction-diffusion master equation (RDME) [6 describes the dynamics of the molec- 
ular copy numbers under the stochastic reaction-diffusion process. In this formulation, the 
spatial domain is partitioned into small voxels so that the molecules are well-mixed inside 
each voxel, but the copy numbers of molecules are small enough to preserve the stochastic- 
ity. Thus, the dynamics of reactions inside each voxel is governed by the chemical master 
equation (CME) [4 , while the diffusion is modeled as the exchange of the molecules between 
voxels. 

In [7 , the stochastic pattern formation was studied for Brusselator [8j based on the RDME. 
Interestingly, it was observed that Brusselator can exhibit spatial patterns even when the 
deterministic reaction-diffusion model converges to a homogeneous equilibrium. This implies 
that intrinsic fluctuation drives a particular spatial mode, and results in inhomogeneous 
spatial structure. This research thereafter inspired the investigation of the noise-induced 
Turing patterns in Levin-Segel model [9 and the noise-induced spatio-temporal oscillations 
in Brusselator [10 . 

In these papers jll [H [10] , an approximate model obtained by linear noise approximation 
(LNA) pT] was used to compute the spatial power spectrum of the patterns. Later, the 
LNA for reaction-diffusion systems was formulated in a more general form in [12^, and 
the analytic results obtained by the approximation agreed with the stochastic simulations. 
Although these methods were shown to give a good approximation, the mechanistic basis of 
the stochastic pattern formation is not necessarily revealed. Hence, it is desirable that the 
system be analyzed from a system theoretic viewpoint to understand the mechanism behind 
the phenomenon. 

The goal of this paper is to provide a mechanistic understanding of the stochastic reaction- 
diffusion system, and to reveal the mechanism of the noise-induced Turing pattern from a 
control theoretic viewpoint. Specifically, it is shown that the computation of the covariance 
of intrinsic noise can be viewed as the I-L2 norm computation of a coupled dynamic multi- 
agent system, where disturbance inputs are injected to each agent's states and sensed out- 
puts. Then, an efficient method to compute the spatial power spectrum is derived by using 
the characteristic structure of the reaction-diffusion system. Using these analytic tools, we 
analyze noise-induced Turing patterns observed in the celebrated Gray-Scott model p!3l[T4] . 
It should be remarked that to the authors' knowledge, this paper gives the first presentation 
of the noise-induced Turing pattern in the Gray-Scott model [131 [14]. Finally, we reveal the 
mechanism of the noise-induced pattern formation with the 7^2 norm interpretation of the 
multi-agent system. 

This paper is organized as follows. In the next section, the RDME is briefly presented. 
Then, an approximated system of the RDME is derived by using the LNA in Section III. 
In Section IV, the characteristic structure of the approximated system is revealed, and its 
control theoretic interpretation is presented. Furthermore, an efficient method to compute 
spatial power spectrum is obtained. In Section V, noise-induced Turing patters are demon- 
strated with the Gray-Scott model, and their mechanism is illustrated. Finally, the paper 
is concluded in Section VI. 

2 Model description of stochastic reaction-diffusion sys- 
tems 

In this section, the dynamics of stochastic reaction-diffusion system, which we call reaction- 
diffusion master equation (RDME) [6 is introduced. 
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Figure 1: The reaction-diffusion scheme considered in the RDME formulation. 



Consider a set of chemical reactions that consists of n molecular species, Mi,M2, - ■ ■ ,Mn, 
and mo reactions, 7?.i, 7?.2, • • • , Ti-mo^ i'^ ^ spatial domain. Suppose the domain is partitioned 
into N voxels, Vi, V2, • • • , Vn with the same volume VL. It is assumed that molecules are 
well-mixed within each voxel, and can react with those in the same voxel. Figure [U shows 
an example of the situation for one dimensional case. 

Let an integer vector Xi G Z" denote the copy number of molecular species inside the 
voxel Vi (i = 1, 2, ■ ■ ■ , N). Stacking the vector Xi, we define 

X:=[Xj,X^,---,XireVf. (1) 

It is known that the stochastic time development of the molecular copy number follows the 
chemical master equation (CME) for well-mixed chemical systems [4^- Thus, the dynamics 
of chemical reactions inside the voxel Vi is given by 

^^^^ = £ {Wr{Xi - Sr)P{Xi -Sr,t) 

r=l 

-Wr{Xi)P{Xi,t)) = AiP{Xi,t), (2) 

where P{Xi^ t) denotes the conditional probability that there are Xi molecules in Vi at time 
t for a given initial state and timeEl- The function Wr{-) : (r = 1, 2, • • • , mo) and 

the vector Sr G Z"^ denote the propensity function and stoichiometry for the reaction TZr 
(see [15 for details). We define Ai as an infinitesimal generator describing the development 
of P(X-,t). 

The diffusion of molecules can be modeled by the exchange of molecules between voxels (see 
Fig. [1]). When one molecule of A4k moves from the voxel Vi to Vj, the number of molecules X 
changes to [Xf , • • • ,Xf -e^, • • • ,Xj + e^, • • • where e/c := [0, • • • , 0, 1, 0, • • • , 0] G 

is the standard unit vector with 1 at the k-th entry. Let li {i = 1, 2, • • • ^N) denote a 
set of adjacent voxels of Vi, i.e., 

•= {j e {1, 2, • • • , A/'} I Vj is adjacent to V^.}. (3) 

In general, X is updated as 

X^Vij^Bk, (4) 

^More precisely, P(Xi,t) should be written as P(Xi,t\XiQ,to) with the initial state Xi^ and time to. 
In this paper, however, we omit the condition, and simply write P(Xi,t) to avoid notational complexity. 
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for each diffusion event, where the vector Vij G has —1 at the i-th entry, +1 at the j-th 
entries (j G li) and at the other entries. Thus, the dynamics of diffusion is written as 

^^^^ = E E E mx.+i)p{x-vi, ® e,,t) 

i=l jeXi k=l 

-dkX,P{X,t))=VP{X,t), (5) 

where t) is the conditional probabihty that there are X molecules at time t for a given 

initial state and time. The diffusion constant of the molecule A1/c is defined by Sk := dk/i'^ 
with a deterministic diffusion rate dk > and the characteristic length of each voxel £. The 
probability that one molecule of Aik moves from Vi to an adjacent voxel within a small 
time interval [t, t + At] is given by 6kXiAt. We define V as the infinitesimal generator that 
accounts for diffusion. 
From ([2]) and (|5j), we have the following reaction-diffusion master equation (RDME). 

^^i^ = J2AiP{Xi,t) + VP{X,t). (6) 

Although the RDME describes the dynamics of intrinsic noise in stochastic reaction- 
diffusion systems, it is known that the RDME is hard to solve analytically in most cases. 
Hence, in the next section, we introduce an approximation method to understand the prop- 
erties of the intrinsic noise in an analytic way. We hereafter restrict our attention to one 
dimensional spatial domain for simplicity. The theoretical results, however, can be extended 
to higher dimensional cases in a similar approach. 



3 Linear noise approximation 

In this section, we briefly describe the idea of linear noise approximation [TT], and derive 
the approximated dynamics of (j6j). 

Let := X/Q G denote the concentration of molecules. Note that the molecular 
copy number X divided by the volume Q is the concentration. It is known that x^ converges 
to a deterministic solution in the thermodynamic limit ft ^ oo [161 Hi] • Specifically, let 

(i = l,2,--- ,iV), (7) 



and X := [xj^ x^ ^ • 
diffusion equation 



Xi 



lim 



Then, x(t) follows the spatially discretized reaction- 



{L®D)x 



(8) 



with the initial value ic(0) := limct- 
matrix 



L := 



-1 1 
1 -2 1 



ic^(O), where L is defined by a discretized Laplacian 




T>NxN 



(9) 
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and D is the deterministic diffusion rate matrix defined by 

D:=dmg[di,d2r-- ,^in] ^R"^" (10) 

with di := Sii'^. In ([9]), the Neumann boundary condition with zero flux is assumed. Though 
the fohowing argument can also be applied to the case where the boundary is periodic, we 
hereafter assume the Neumann boundary for simplicity. 

As we have seen in ([7]), the variance of molecular copy number Xi becomes less significant 
as the volume Q becomes larger, and eventually converges to zero. Thus, one would expect 
from this observation that Xi is approximated around the deterministic solution Xi as 

Xi Qxi + J^^i (i = 1, 2, • • • , TV), (11) 

where r]i is a noise term whose properties are specified later. It should be noted that Qxi is 
dimensionless, and the noise grows with 0(1^^/^). 
Let 77 e M^^ be defined by 

rj:=[rjl,rjl---,r^%feR-'^. (12) 

Note that the subscript of r] stands for the index of a voxel. Since Xi is a determinis- 
tic value obtained from (|8]), 77 is the only random variable that determines the stochastic 
fluctuations. Thus, defining n(?7,t) as the probability distribution of 77 at time t, we have 
P(X,t) r^^~n(r7,t), where the scaling term comes from the change of variable of prob- 
ability distribution. We see that the probability distribution 11(77, t) contains the property 
of intrinsic noise. 

The idea of linear noise approximation [TT] is that we replace X in the equation (|6]) with 
the right-hand side of (pT]) . and expand around the deterministic value Qx with Taylor 
expansion. Sorting with the power of we have the deterministic model (|8j) in the term 
of r^^/^, and a Fokker-Planck equation of 11(77, t) ^he term of ft^. Since 0{ft~^^'^) terms 
become less significant as Q becomes large, we can approximately adopt the Fokker-Planck 
equation as the dynamics of 11(77, t). Then, the standard argument of stochastic systems 
allows us to see that 77 satisfies the following dynamic equation. 

dv = Mt)ridt + Sw'^^^^dB(t), (13) 

where B(t) is a vector oi R \= moN + min independent standard Wiener process with 
mi := J2f=i The matrix J^^t) ^ R^^^^^ is Jacobian of (|8|) at x{t). S and W^^t) 
are defined as stoichiometry matrix and a diagonal matrix with the propensity functions at 
diagonal entries, which are specified from ([6]), respectively [15]. 

In the following section, we consider the stochastic fluctuation around a spatially homo- 
geneous equilibrium point, at which stochastic Turing pattern is expected. In particular, a 
characteristic structure of the stochastic system (fT3)) is revealed based on properties of the 
reaction-diffusion system. 

4 Analysis of noise-induced spatial patterns 

4.1 Structure of the stochastic system and noise covariance 

Let Xe G denote an equilibrium of Xi = fixi)^ where /(•) is given in the deterministic 
model defined in ([8|). We can see that x* := [x^^x^^--- ^x^]^ G is a spatially 

homogeneous equilibrium point, i.e., the values are the same between voxels. 
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Table 1: Meaning of the matrices 



5*0 Stoichiometry matrix for T^i , 7^2 , • ' ' '^rj 



Si Incidence matrix of voxels 



Wq Diagonal matrix whose entries are propensity functions for 
Wi Diagonal matrix whose entries are the equilibrium concentrations of 
D Diagonal matrix whose entries are the diffusion rate of A/i , A/2 , • • • , Mn • 
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Figure 2: Block diagram of the stochastic system (p!3|) . 



Using the structure of (|6j), it can be shown that the matrices in ([13]) are given by 



Jx* := In ® K ^ L D e 



r)nNxnN 



S:= [In® So, Si ^ In] G Z"^><^, 
In ®Wo O 



O 



In 



DWi 



)RxR 



where 



So 
Si 
Wo 
Wi 



dl 
dx 

[«Sl,«S2, 
[Vl,V2, 

= did.g{wi{Xe),W2{Xe)r ■ ■ ,^mo(a?e)) ^ 

= diag(x(i\x(2\--- G R^><^. 



x=x* 



The vector Vi is := [vij^^Vij^^--- ^Vij^^,^] G Z^^l"^^' with jk G X^. The function Wi{-) : 

R!f: ^ is defined as the deterministic reaction rate for reaction and Xe*^ is the i-th 
entry of cCg (i = 1, 2, • • • , n). The meanings of each matrix are summarized in Table [TJ 

Substituting the above definition into (p!3|) , we see that the overall system has the structure 
shown in Fig. O Since the matrices ^o, Wo and K are associated with the reactions inside 
each voxel, Bq excites the intrinsic noise arising from intravoxel reactions. On the other 
hand, Bi excites the one arising from diffusion (see Fig. [2]). 
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It follows from (p!3|) that the covariance of the noise E[r7?7^] the steady state cc* is given 
by the Lyapunov equation 

J^*S + SJJ. + SW:,^S^ = 0. (14) 
Thus, we have the following proposition. 

Proposition 1. Consider the system [T^} . Suppose the spatially homogeneous equilibrium 
X* is locally stable. The steady state covariance of the noise S := E[?7?7-^] is given by the 
positive definite solution of the Lyapunov equation 

(/at + L D)^ + ^{In ^L^D) 

+ (/Ar ^ SoWoS^ DWi) = (15) 



This proposition provides an approximation of the covariance of the intrinsic noise whose 
dynamics is given by ([13]). The i-th n by n diagonal block of S stands for the covariance 
inside the voxel Vi, i.e., ¥,[r]ir]f]. The covariance between voxels appears in off-diagonal 
entries, where the {i^j) off-diagonal block is defined by ¥,[rjir]j]. It should be noted that 
Tr(S) in Proposition 1 provides I-L2 norm of the system in Fig. [2l 

Remark 1. It is interesting to observe the system in Fig. [2] from a viewpoint of dy- 
namical multi-agent systems. In Fig. [2l the blocks in the reaction part (upper blocks) are 
homogeneous block diagonal matrices, and each block diagonal entry H{s) := {si — K)~^ 
corresponds to the linearized dynamics of intravoxel reactions. It is clear that Bo is the 
noise that perturbs the states of each subsystem H{s). The blocks in the diffusion part 
(lower blocks), on the other hand, describes the structure of information exchange between 
H{s). Since is the incidence matrix associated with the graph Laplacian L, Bi can be 
considered as a perturbation to the sensed output ry^-i— 77^. Thus, one might expect that 
the study of this class of system is useful for engineering applications as well. □ 

4.2 Spatial power spectrum analysis of the intrinsic noise 

Although the matrix S in ([15]) contains all information about the spatial covariance of 77, it 
is not easy to see the existence and profiles of spatial pattern at a glance. Hence, we here 
consider spatial spectrum analysis based on Proposition 1. 

The following theorem presents an efficient computation method of the spatial power 
spectrum of the noise. 

Theorem 1. Consider the system / fl^) . Suppose the spatially homogeneous equilibrium 
X* is locally stable. Let (z = 1, 2, • • • , A/") denote the spatial frequency components of the 
steady state noise rji, i.e. 




Then, the power spectral density := E[^/c^^] G W^^'^ at the frequency uok •= (^ — 
l)7r/A/' (/c = 1, 2, • • • , N) is given by the positive definite solution of 

-2XkDWi = 0, (17) 
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Figure 3: The decomposed subsystem obtained from the large system in Fig. O 



where 

A.:=-4sin^(fcil^). (18) 

We see from Theorem 1 that the spatial power spectral density of the molecule M.i is 
obtained as the (z,i)-th entry of the matrices (/c = 1, 2, • • • ^N). In fact, S is obtained 
by applying discrete cosine transform to the covariance matrix S. It should be noted that 
S/e (/c = 1, 2, • • • , A^) are obtained by solving N Lyapunov equations of the size n by n, 
which is more computationally efficient than solving nN by nN Lyapunov equation ([T5]) 
and applying discrete cosine transform. 

From a control theoretic viewpoint, the transformation presented in Theorem 1 cor- 
responds to the decomposition of the system into N subsystems as shown in Fig. [3l 
where A = Ai, A2, • • • , Aat. We see that the feedback gain A corresponds to frequency 
uok {k = 1,2, •• • , A/"), and the spatial spectral density varies in terms of A. Thus, the 
computation of the spatial power spectrum is essentially the same as I-L2 norm computation 
for various feedback gains A. 

Since the noise Bo and Bi are decoupled in Fig. [H the I-L2 norm can be independently 
computed for each input. Thus, we have the following proposition. 

Proposition 2. The Lyapunov solution Ek in ([7?| j can be decomposed into Ek = S^i + 
S/c2 = I7 2, • • • , A^); where Ej^i and are the positive definite solutions of 

{K + XkD)Eki + {K^ + XkD) + SoWoS^ = 0, (19) 
{K + XkD)Ek2 + Ek2 {K^ + XkD) - 2XkDWi = 0. (20) 



It should be noted that Eki and Ek2 represent the spectral density of the noise originating 
from intravoxel reactions and diffusion, respectively. Thus, Proposition 2 implies that we 
can independently analyze the contribution of these noises. 

Remark 2. Theorem 1 and Proposition 2 can be applied for systems with the periodic 
boundary condition. In that case, the Laplacian eigenvalues should be alternatively defined 
as Xk = —4 sin^((/c— l)7r/A/'), and the solution E corresponds to the discrete Fourier transform 
of S. □ 



8 



Table 2: Reactions of the Gray- Scott model 





Reaction 


Propensity Wi{-) 
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5 Numerical Simulations 

In this section, we first illustrate stochastic Turing patterns induced by intrinsic noise with 
the well-known Gray-Scott model [131 [14]. Then, we show that the spatial power spectrum 
analysis presented in the previous section can capture the noisy spatial patterns. Fur- 
thermore, we present underlying mechanism of noise-induced Turing patterns based on the 
system theoretic interpretation given in the previous sections. Due to the limitation of the 
space, some details are left for our future publication. 



5.1 Noise-induced Turing patterns with the Gray-Scott model 

The Gray-Scott model consists of the five chemical reactions in Table [21 In reaction T^i, 
U is converted into where V catalyzes the production of itself. Then, the reaction IZ2 
changes V into the final product form P. The reactions 1Z2,^1Za and IZ^ describe the inflow 
of U from outside, degradation of U and degradation of V. 

The deterministic reaction-diffusion model of the Gray-Scott model is given by 

^ = —uv'^ + a(l — u) ^ dV'^u 

f (21) 

— = uv^ - (a + k)v + V V, 

where u G M+ and v G M+ are normalized concentrations of U and V . The dimensionless 
constants a, k and d are normalized reaction rates and a diffusion rate defined with the 
reaction constants in Table [2 as 

KiUq kiUQ d2 

Let the parameters be set as ('Uo, ^1, ^2, ^a, <^2) = (3.0,4.0 x 10~^, 1.98 x 10~^,2.16 x 
10-^ 2.16 X 10-^^ 3.6 X 10-^"^), which corresponds to (a, 6, d) = (6.0 x 10-^ 5.5 x 10-^ 6.0). 
With this parameter set, both deterministic and stochastic simulations exhibit the spatial 
pattern as shown in Fig. [H For the stochastic simulation, the Gillespie's algorithm [18] was 
used, where the domain was partitioned into = 32 voxels, and the characteristic length 
and the volume were I = 1.0 x 10~^ and Q = 1.0 x 10^, respectively. 

We now change the parameter /c2 in the above example to 1.76 x 10~^, which results 
in (a, 6, = (6.0 x 10~^,4.9 x 10~^,6.0). According to the deterministic analysis [19 , a 
spatially homogeneous equilibrium is stable, hence no spatial pattern is expected. In fact, 
the deterministic model converges to a homogeneous equilibrium (Fig. [5l(Left)). However, 
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Figure 4: Time development of the concentrations/copy numbers of U for the first parameter 
set. (Left) A simulation result of the deterministic model. (Right) A result of the stochastic 
simulation. 




Figure 5: Time development of the concentrations/copy numbers of U for the second pa- 
rameter set. (Left) A simulation result of the deterministic model. (Right) A result of the 
stochastic simulation. The red and blue patches imply the existence of noise-induced spatial 
pattern. 



we observe that the stochastic simulation in Fig. [5] (Right) still displays the spatial pattern. 
It should be emphasized that the deterministic model (|2T]) failed to capture this spatial 
pattern. 

5.2 Spatial power spectrum analysis 

In this section, we first show that Theorem 1 and Proposition 2 can capture the stochastic 
spatial patterns shown in Fig. [5] (Right), then illustrate the mechanism of the stochastic 
pattern generation from a system theoretic point of view. 

Consider the case where the deterministic system (|2T]) does not present Turing pattern. 
Using Theorem 1, we can approximately compute the steady state spatial power spectrum. 
Figure [6] shows the computation result in terms of b. We see that as b becomes larger, 
sharp peak appears at /c = 6 and 7, which correspond to ujq = 57r/32 and ur = 67r/32. For 
5 = 4.9x10"^, which is the second parameter set in the previous example, the peak frequency 



10 




10 20 30 

Frequency {k) 



Figure 6: Spatial power spectrum of intrinsic noise rj computed from Theorem 1. The 
horizontal axis is k m ujk = — l)7r/A^. 
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Figure 7: Trajectory of the eigenvalues oi K ^ XD in terms of A. The mark * stands for the 
eigenvalues at A = A/^. 

is UJ7. Hence, we can expect spatial patterns with three periods. Indeed, the spatial structure 
can be confirmed in Fig. [5] (Right). The ingredient of the power spectrum can be analyzed 
with Proposition 2, and it can be concluded that both reactions and diffusion contribute to 
form the peak at non-zero frequency in Fig. [6] 

It is interesting to consider these observations from a control theoretic viewpoint. In 
general, the deterministic Turing instability is determined from the stability of K -\- XD, 
which corresponds to the stability of the system in Fig. [3l Figure [71 illustrates the eigenvalue 
distribution of K -\- XD in terms of A. Note that this corresponds to drawing the poles of the 
system in Fig. [3] by changing the feedback gain A. We see that one eigenvalue approaches 
to the imaginary axis as k increases, and the largest real part is achieved at /c = 8 (see Fig. 
[7]). In view of the stability of this system, it is natural that the largest I-L2 norm is obtained 
around k = 8. This observation corroborates the aforementioned analysis illustrated in Fig. 
ini where we actually have the peak at k = 7. 
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6 Conclusion 



We have analyzed noise-induced pattern formation under the reaction-diffusion kinetics de- 
scribed by the RDME. Using the hnear noise approximation, we have revealed the charac- 
teristic structure of the stochastic reaction-diffusion system, and presented the underlying 
mathematical structure of noise-induced pattern formation from a control theoretic point of 
view. The analytic tool to examine the spatial power spectrum of intrinsic noise has also 
been derived in the analysis, and its effectiveness has been confirmed through the numerical 
simulations. 
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